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Abstract. Two algorithms that combine Brownian dynamics (BD) simulations with mean-field 
partial differential equations (PDEs) are presented. This PDE-assistcd Brownian dynamics (PBD) 
methodology provides exact particle tracking data in parts of the domain, whilst making use of 
a mean-field reaction-diffusion PDE description elsewhere. The first PBD algorithm couples BD 
simulations with PDEs by randomly creating new particles close to the interface which partitions 
the domain and by reincorporating particles into the continuum PDE-description when they cross 
the interface. The second PBD algorithm introduces an overlap region, where both descriptions 
exist in parallel. It is shown that to accurately compute variances using the PBD simulation requires 
the overlap region. Advantages of both PBD approaches are discussed and illustrative numerical 
examples are presented. 
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1. Introduction. Spatial reaction-diffusion models have been widely used for 
the description of biological systems [25] . Often continuum approaches, written in 
the form of reaction-diffusion partial differential equations (PDEs), are used due to 
their simplicity and the vast number of ready-to-use numerical solvers. However, 
many biological effects cannot be fully described by deterministic PDE-based models. 
This is because a deterministic model requires large copy numbers of molecules to 
minimize the relative fluctuation of the spatial concentration. If low copy numbers 
are present in a biological system [55] , then stochastic models such as mesoscopic 
compartment-based algorithms [191 [S] or trajectory tracking (Brownian dynamics) 
methods, may be deployed [2| [29] . 

In many situations individual trajectories are important only in certain parts of 
the domain, whilst in the remainder of the domain a coarser, less detailed, method can 
be used [T5]- This is the case, for example, in the modelling of ion-channels [M]- Ions 
pass through a channel in single file and an individual-based model has to be used to 
accurately compute the discrete, stochastic, current in the channel [J. The positions 
of individual ions are less important away from the channel where copy numbers 
may be very large (rendering a detailed Brownian dynamics description infeasible) 
[5] . Another example is the stochastic reaction-diffusion modelling of filopodia which 
are dynamic finger-like protrusions used by eukaryotic motile cells to probe their 
environment and help guide cell motility [33] . These relatively small protrusions are 
connected to a larger cytosol compartment. If a modeller is interested to understand 
the dynamics of filopodia, then there is a potential to decrease the computational 
cost of simulations by using a coarser model in the cytosol. In both examples, it is 
important to understand how models with a different level of detail can be used in 
different parts of the computational domain 15J. 

In this paper, we develop algorithms that calculate Brownian dynamics (BD) 
paths in a desired part of the domain, whilst using a continuum PDE-based model 
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in the remainder. This PDE-assisted Brownian dynamics (PBD) methodology has 
the advantage that efficient methods for solving PDEs can be used for large parts of 
the modelled domain, whilst BD data is available in other areas where required. The 
main goal of the PBD methodology is to get the same statistics (means and variances) 
in the BD subdomain as we would get if we we were able to use BD simulations in 
the whole domain. In particular, the correct coupling between the two parts of the 
domain is of vital importance for the accuracy of a PBD algorithm. 

The paper is organised as follows. Section [2] states, in mathematical terms, the 
requirements for the developed algorithms and introduces the notation used through- 
out the remainder of the paper. We then introduce the first PBD algorithm for a 
pure diffusion system in Section [3j where we also explore the complications of this 
algorithm. In Section [4] we present the second PBD algorithm which provides more 
accurate computations. In Section [5] we investigate issues relating to the introduction 
of reactions into the system and present several computational examples. 

2. Problem formulation. Consider a general Brownian dynamics reaction- 
diffusion simulation with M chemical species in the (open) domain J] C t 3 . We 
denote by ?ij(x, t), j = 1, 2, . . . , M, the expected spatio-temporal concentration of the 
j-th chemical species at the position x and time t over our domain SI. The approxi- 
mate mean-field reaction-diffusion PDEs for the time evolution of concentrations can 
be written as follows 

^=D j Ap j + R j (p 1 ,p 2 ,...,p M ), j = l,2,...,M, (2.1) 

where pj = pj(x, t) : SI x [0, oo) — » [0, oo) is the mean-field approximation of rij, Dj is 
the diffusion constant of the j-th chemical species and Rj : [0,oo) M — > K represents 
the reaction terms. 



The goal of PBD algorithms is to couple macroscopic description ( 2.1 ) in the open 
subdomain Sip C SI with a stochastic BD simulation in the open subdomain Sip C SI, 
where the closures of fip and Sip cover SI, i.e. 

Slcn^USTp". (2.2) 

In Sip, we will consider BD trajectories of individual molecules, i.e. the state of the 
microscopic subdomain Sip is defined by the number Ng\t) of molecules of the j- 
th chemical species at time t and their positions xf\t) € Sip, i = 1,2,... ,Ng\t), 
j = 1,2,..., M. We denote by I the interface between the subdomains Sip and Sip, 
namely 

I = Ti~B~n TLp~. (2.3) 

In this paper, we will investigate two cases: 

[A] fls and Sip do not overlap, i.e. Op n J7p = 0; 

[B] there exists an overlap region where the PDE description and BD simulations 
exist in parallel, i.e. Sip n Sip ^ 0. 



The case [A] will lead to the PBD algorithm (Al)-(A5) presented in Table 3.1 The 



case [B] is implemented in the second PBD algorithm (B1)-(B5) which is presented 



in Table 4.1 We will start our discussion with the case [A] because it is less technical 



to implement than the case [B] 
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Fig. 2.1. Sketch of the first PBD algorithm and the notation related to it. In Qp, molecules 
are described by their density distribution p(x,t) and in the microscopic domain Qb described by 
the number iVg(t) of molecules and their positions Xi(t), i = 1, 2, . . . , Ng(t). The interface between 
these domains is denoted I. 



To simplify our presentation, we will consider that Q is a "narrow" three-dimen- 
sional domain, and hence only consider the process mapped onto an effective one- 
dimensional domain 57 C R by assuming that the system is well mixed in the other 
two dimensions. In particular, we have fla C C 1 and the state of the BD 
sub-domain fl B will be described by the x-coordinates of molecules which we will 
denote as x\ j \t) G Q, B , i = 1,2,..., N^\t), j = 1,2,..., M. We simulate the system 
using finite time step At > 0, in which particles in Q B change their position according 
to the discretized version of the overdamped Langevin equation 



,r 



U \t + At) =x<f\t) + ^2D^Kti, (2.4) 



where £ is a normally distributed random variable with zero mean and unit variance. 
Figure [2~T] shows a sketch of the described system for one chemical species along with 
the notation used in the case [A] . Our goal is to construct PBD algorithms which will 
satisfy the following two conditions: 

Condition (C.l): We require that the expected distribution of molecules in il B \H,p 
match that of the expected distribution rij(x,t), i.e. the distribution which we would 
obtain if we used detailed BD simulations in the whole domain fi. In particular this 
means that for every set A C £l B \ Qp, the expected number of particles in A at time 
t > has to satisfy 

E {^(i) G A,i = l,...,N i J ) (t)} =J nj (x,t)dx, \/Acn B \n P . (2.5) 

The choice of an arbitrarily small but finite interval A = [x,x + dx) for x G £Ib \ 
leads to an alternative formulation of this condition 



E 



(x^ (t) G [x,x + dx) ,i = 1, . . . , Njp (t)\ =n j (x,t)dx, yxefl B \n P . 

(2.6) 

Condition (C.2): Whilst we aim to match the expected outcome of the stochastic 
simulation to n(x, t), we also want the variances of molecule distribution in £l B \ 
to match that which would be expected if a BD simulation were to be performed over 
the entire domain fl. 

In fip \ £l B , the system is described by the concentration vector pj, j = 1,2,..., M, 



which evolves according to the PDE (2.1). This distribution, whilst continuous, is 
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not strictly deterministic since it is coupled with the stochastic outcomes of the BD 
subdomain Sip. If the stochastic reaction-diffusion model only includes zero-order 



or first-order reactions, then the mean-field PDE (2.1 1 describes the expected be- 
haviour of stochastic models [T^] . In this case it is reasonable to require the following 
additional condition. 

Condition (C.3): We require 

e [pj(x, t)] = n,j(x, t) , y x g si P \ sip , t > o. (2.7) 



In the case [A], the conditions (C.1)-(C.3) can be simplified by observing that SI b = 
Sip \ Sip and Sip = Sip \ Sip. In the case [B], we will also require that the PBD 
algorithm gives the correct mean distribution of molecules in the overlap region O = 
Sip n Sip, i.e. 

E\pj(x,t)]dx+E (x\ j) (t) G [x,x + dx) ,i = l,...,^'^)} =nj(x,t)dx (2.8) 

for all x G O = Sip n S1 P and t > 0. 

3. PBD simulation of diffusion. In this section, we explain our case [A] PBD 
algorithm (i.e. SlpflSlp = 0) using a system of diffusing non-interacting molecules of a 
single chemical species. Therefore, for all further discussions in this and the following 
section we will drop the index j representing the species. Then the macroscopic PDE 



(2.1) for this system becomes the diffusion equation 



where n = n(x,t) : SI x [0,oo) — > [0,oo) and D is the diffusion constant. We will 
consider the infinite domain SI = R for simplicity. Without loss of generality, we 
assume that Sip = (— oo,0) and Sip = (0,oo), i.e. the internal boundary / = {0} 
is situated at the origin x = 0. In the case of diffusion only, the total number of 
molecules N in the system is conserved, i.e. 

N= / n(x,t)dx, foralU>0. (3.2) 
Jn 

Hence, for the PBD algorithm the conservation of mass condition takes the form 

N= [ p{x,t)dx + N B (t) , for alii >0. (3.3) 

Since the diffusing molecules are non-interacting, we can express Condition (C.2) in 
mathematical terms. If all particles start with the same initial condition, then each 
particle has the (identical) probability — f A n(x,t)dx/N of being in a set A at 
time t. Consequently, the expected number of particles in A at time t is NpA and the 



variance is equal to Npa(1 — Pa)- Substituting = f A n(x,t)dx/N and using (3.2), 
we get Condition (C.2) in the following form 

vax[\{ Xi (t)eA,i=l,...,N B (t)}\]= f n(x,t)dx(l~^p^) , (3.4) 

J A V J n n{x,t)dxJ 



■5 



p(x,t) 
p(x, t + At) 




Q P I 


iCt(t + At) 



X 

Fig. 3.1. Illustration of step (Al) of the first PBD algorithm. Dashed line: p(x,t); solid line: 
p(x, t + At) as defined in \3.6\; shaded area: ct(t + At) as defined in \3.%\ . 



for all A C Sip and t > 0. Using again A = [x,x + dx), we obtain the alternative 
formulation 

var [\{xi(t) e [x,x + dx) ,i = 1, . . .,N B {t)}\] = n(x,t)dx . (3.5) 



In Sections |3.1| and |3.2| we present one update step of the continuum and the 
particle-based simulations respectively, before the full PBD algorithm (A1)-(A5) is 
formulated in Section [3~3| In Section [3~4} we will discuss the accuracy of the algorithm 
with respect to the Conditions (C.1)-(C.3). 

3.1. Updating the PDE regime in Sip. At time t, we have the concentration 
p{x,t) for x G fip and are aiming to calculate the concentration p(x,t + At) that 
corresponds to a realisation of one time step At of the diffusion process ( |3.1[ ). We 
therefore define the exact outcome of a diffusion step in the full domain SI given initial 
data p(x, t): 

p(x,t + At)= K(x-x',At)p(x',t)dx' , (3.6) 
J n P 

where K(x — x' , At) is the diffusion kernel 

^■"'-Tara-" ("*&)• (3J) 

and the function p{x, t) has support Q. Using this process, a certain proportion of 
the concentration distribution, namely 

a(t + At)= / p{x,t + At)dx (3.8) 



would have crossed the interface / in the time interval [t, t + At) (see Figure 3.1 1. This 
value a(t + At) represents the expected number of molecules to cross the interface 
from Hp to VIb in the time interval [t,t + At). Since all molecules in Sip are identical 
(of the same chemical species and have the same probability distribution proportional 
to p(x,t + At)), the number of molecules that cross the interface from Sip to Sip in 
the time interval [t,t + At) is Poisson distributed with average a(t + At). Let us 
assume that At has been chosen small enough to ensure that a(t + At) -^i 1. In this 
case the probability of more than one particle crossing the interface is negligible and 
we need to consider two cases: 



() 



(i) one particle gets created in Op with the probability a(t + At); 

(ii) no particle is created with the probability 1 — a(t + At). 

For both cases we need to calculate the updated concentration p(x, t + At) where 
x € Op. 

We consider the concentration p(x,t) as the distribution of N — Ap(i) identically 
distributed particles at time t. Therefore each of these particles has at time t the 
probability distribution p(x, t) / (N—Nb (t) ) . After one time step each of these particles 
can be found in an infinitesimal interval [x, x + dec) with probability p(x, t + At) dx. 
For each particle, its probability distribution given that it did not leave Op can be 
calculated as 

, p(x,t + At) 

t + M) = N-N B (t)-a(t + At) ' f ° r X G QP • 

On the other hand, if a particle does leave the domain Op, then its distribution 
function becomes zero for x £ Op. Instead the particle is introduced into Op at time 
t + At at a location x with the probability distribution given by 

p 2 (x,t + At) = P ( X ;* + *f ) , forzeOp. (3.9) 
a(t + At) 

This updating process is rather like collapsing a wavefunction in quantum mechanics 
[30] • We have a look to see if a particle has crossed the boundary into Op: if it is 
there on the other side then its distribution function collapses to a S function at its 
new position, while if it is not, then the distribution function collapses to zero in Op 
with corresponding rescaling in Op. 

Combining these arguments for all the particles, we see that the last update step 
is a simple rescaling of the probability distribution so that the updated distribution 



satisfies conservation of mass according to (3.3). We therefore have 



/ » ,\ f Pa) p(x,t + At), in the case (i), . „ ,„ 
t + ^ = { '^x, t + At), in the case (ii). f ° r X S ' t 3 ' 10 ) 

with the rescaling constants j3u) , given by 

= N-N B (t)-i(t + At) ' = N — N B (t) — a(t + At) " (3 - H) 

Note that the update step for the continuum regime satisfies conservation of mass 



(3.3). An illustration of the two cases can be seen in Figure 3.2 



3.2. Updating the BD regime in Op. We use the discretized version of the 



overdamped Langevin equation introduced in (2.4) to update the positions of particles 
in Op. If the position of the i-th molecule, computed by (2.4), is in Op at the end 
of the time step, then we will continue representing it as a particle. Note, that a 
particle that crossed the boundary I and came back into Op during the time step 
[t, t + At) is also captured by this case. We have to be more careful whenever the 



position Xi(t + At), computed by (2.4), is inside the PDE subdomain Op at time 
t + At, i.e. Xi(t + At) e Op. In this case, the random walk of the i-th molecule 
crossed the interface I without crossing back at the end of the time step. This event 
needs to be taken into account for the update of the final concentration p(x, t + At) 
in Op. As we know the exact position of this particle i at time t + At we add a Dirac 
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p{x, t + At) j N B (t) + 1 p(x,t + At) j N B (t) 
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n, 



X X 
(a) Case (i) (b) Case (ii) 

Fig. 3.2. Illustration of steps (A2) and (A3) of the first PBD algorithm. Dashed line: p(x,t- 
At); solid line: p cont (x,t + At); (green) circle: created particle. 



5 function at the position Xi(t + At) S rip. Therefore we compute p{x, t + At) in fip 

by 

p(x,t + At) =p cont (x,t + At) + S(x-Xi(t + At)), 



where p CO nt(x,t + At) is given by (3.10). 



3.3. The first PBD algorithm. This algorithm computes the concentration 
p(x,t) for x € rip, and the number N B (t) and positions of BD particles Xi(t) € flp, 



i = l,2,..., N B (t). One time step of the first PBD algorithm is presented in Table 3.1 
as algorithm (A1)-(A5). In order to simplify the presentation of this algorithm, we 
consider that the time step At is chosen so small that a(t + At) <C 1. In particular, we 
only need to implement cases (i)-(ii) presented in Section [3. 1| because the probability 
that two or more molecules are initiated in 51 p during one time step is negligible. 

The auxiliary distribution p(x, t + At) in step (Al) is in practice calculated using 
a numerical approximation algorithm. To calculate a(t + At) we can either use this 
numerical approximation of p[x, t + At), which requires an additional time step At <C 
At to be used to ensure the accuracy of ait + At), or (more efficiently) we can 
approximate a(t + At) analytically using a boundary layer expansion in the vicinity 
of the interface /. 

By construction, the algorithm (A1)-(A5) satisfies the conservation of mass con- 



dition (3.3 1 and the concentration p(x,t) satisfies non-negativity. We will now show 
that this algorithm also guarantees the correct expected outcome and therefore satis- 
fies Conditions (C.l) and (C.3). 

Theorem 3.1. Consider the BD simulation of N diffusing molecules in the 
computational domain Q which is divided into sudomains 51p C 51 and Qp C fi 
satisfying ( |2.2[ ) and the case [A]. Suppose that N B {0) particles are initially in Qb 
at positions Xi(0), i = 1, 2, . . . , N B (0). Let us initialize p(x,0) as sums of Dirac 5 
functions describing molecules which are initially in Hp, i.e. p(x,0) — n(x, 0) for 
x G fip. Then the expected outcome of the PBD algorithm (A1)-(A5) (presented in 
Table [3T|) satisfies the Conditions (C.l) and (C.3) for arbitrary At > 0. 



Proof. We will show that the identities (12. 51) and (2.7) hold during one iteration 



(A1)-(A5) presented in Table 3.1 It will then follow by induction that they hold for 
all times kAt, k = 0, 1,2, . . . . 

Let us assume that the probability distribution of particles in ilp at time t is 
given through n(x,t), x € flp. Given the probability distribution p{x,t), x € Hp, at 
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(Al) Calculate p(x,t + At) using and a(t + At) using ( |3~8 ). 
(A2) Generate uniformly distributed random number r in (0, 1). 

(i) If r < a(t + At), then create new particle in ftp according to the 
probability density p%(x, t + At) defined in ( |3.9[ ). Set /3 = /3(j) where 
@u\ is given by (3.11). Set Nb,i = 1. 



(ii) If r > a(t + At), then set /3 = /3ui) where (3^ is given by (3.11 ). 
Set N B .i = 0. 

(A3) Compute positions Xi(t + At), i = 1,2, .. . ,Ns(t), of BD particles ac- 



cording to ( 2.4 ). 



(A4) Compute new concentration in ftp by 

p(x,t + At)=/3p{x,t + At) + S ( 

x — Xi(t + At)) , for x £ Ftp . 

i,(i+At)6fi P 

(A5) Update the number of BD particles by 

N B (t + At) = N B (t) + N B>1 - \{ Xi (t + At) £ ftp ,i = l,...,N B (t)}\ . 

Terminate computation of trajectories of BD molecules which landed in 
ftp (i.e. the BD particles which satisfy Xi(t + At) £ ftp.) 
Then continue with step (Al) for time t + At. 



Table 3.1 

One time step of the first PBD algorithm for a system of diffusing molecules. 



time t, the conditional expected value of p(x, t + At) for a; £ ftp is given by 

E [p(x, t + At)\ p(x, t) ]=a(t + At) /%) p(x, t + At) 

+ (1 - a(t + At)) m p{x, t + At) 

+ / K(x- x',At)n(x',t)dx' , 



where and are given by (3.11 1, K(x — x' , At) is the diffusion kernel given in 



( |3.7[) an d the last term represents particles that moved across the interface /, defined 
by ( |2.3[ ), during the time step [t,t + At). Using ( |3.11 1, we obtain 



E [p(x, t + At)\ p(x, t) ] = p(x, t + At)+ K{ x — x , At) n(x' , t) dx 1 . 



Using (3.6 1 and the law of total expectation (law of iterated expectations), we get 

E [p(x, t + At) ] = / K{x-x',At)E[p{x',t)]dx' + K(x - x',At)n(x',t)dx'. 

Using the induction assumption that E[p(ir,i)] = n(x,t), we obtain 

E [p(x, t + At) ] = / n(x',t)K(x - x' , At)dx' = n(x,t + At) , for x £ ftp, 
isi 



i.e. we have derived (2.7) 
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Let us consider a set A C ilp. Given the probability distribution p(x, i), x € fip, 
at time i, the conditional expected number of particles in A at time t + At is 



E [ \{xi(t) e A ,i = 1, . . . , N B (t)}\ | p(x, t) ] = / + At)p 2 (:M + At) da; 

A7a; — x , At)n(a;',t)da;' da;, 



i A Jn B 

where the first term represents newly created particles from the PDE regime fip and 



the second term represents the movement of particles inside f2p. Using (3.6), (3.9) 
and the law of total expectation, we obtain 

E [ | {xi(t) G A , i = 1, . . . , N B (t)} | ] = / / K(x - x', At) n(x', t) dx' dx 

J a Jn 



A 



n(x, t + At) da; , 



which is the condition (2.5). Thus we have showed that both Conditions (C.l) and 



(C.3) are satisfied. This concludes the proof. □ 



Theorem 3.1 also holds if the algorithm (Al)-(A5) is extended to the creation of 
more than one new particle per time step as long as the expected value of the number 
of created particles is a(t + At) and the rescaling is done accordingly. The algorithm 
(A1)-(A5) and the proof can be easily extended for a finite domain fl — [0, L] with 
no flux boundary conditions by redefining the kernel At) accordingly. 



3.4. Discussion of the PBD algorithm (A1)-(A5). In Theorem |3.1| we 
showed that the PBD algorithm (Al)-(A5) satisfies the Conditions (C.l) and (C.3). 
However, we still need to check whether the Condition (C.2) on the variances is also 
satisfied. 

To investigate the variances created by this algorithm, we show the outcome of 
an illustrative numerical example in Figure |3.3| We simulate the diffusion of 100 
molecules in the domain ft = [—1,1] which are initialized at the same location x — 
—0.95. We use no flux boundary conditions. We test the algorithm (A1)-(A5) where 
ftp = (—1,0), £Ib — (0,1) and / = {0}. To calculate p(x, At) we use an implicit 
Euler-scheme with Aa; = 0.01 and a numerical time-step of At = 10~ 6 . The time 
step At used by the algorithm (A1)-(A5) is At = 10~ 3 and we simulate the system 
until Tfinai = 0.2. The time step At <C At for the approximation of p(x,t + At) 
was chosen in order to minimise numerical artefacts. We run 1000 realisations of 
this process and measured the number of particles in 10 intervals ('bins') of the size 
0.1 in Hb- Averaging over 1000 realisations, we calculate the mean value and the 
variance of the particle number for each of the bins at time Tfj na i. The results are 
presented in Figure |3.3| as gray histograms. In this example, it is easy to calculate 
the correct distribution n(x, t) which the algorithm (A1)-(A5) tries to approximate. 



It can be obtained by solving the diffusion equation (3.1) in the domain fl = ( — 1, 1) 



with n(x, 0) = 100 <5(a; + 0.95) and no flux boundary conditions. The expected values 



for both means and variances are plotted as (red) dashed lines in Figure 3.3 

In Figure |3.3| (a) , we see that the mean value of the simulation results matches 
well with the solution of the diffusion process, with only small fluctuations close to 
the internal boundary J at x = due to stochastic effects and inaccuracies caused by 
the numerical approximation of p{x,t + At). However, in Figure 3.3(b)| it becomes 
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(a) Mean (b) Variance 



Fig. 3.3. Simulation results of diffusion of 100 particles in (—1,1) with no flux boundary 
conditions initialized at x = —0.95, i.e. n(x,0) = 100 &{x + 0.95). Results averaged over 1000 
realisations. Dashed (red) line: expected outcome, (a) Solid line: mean value in Qp; (gray) bars: 
particle concentrations in Qb- (b) (Gray) bars: measured variances in particle concentrations. 
Parameter values are described in the text. 



clear that the variance between different realisations is higher than the desired value, 
in particular close to the internal interface. We can explain this effect clearly using a 
thought experiment. 

Let us consider a situation where p(x, 0) is close to the internal boundary / 
and has a peak of mass 99 arbitrarily far away from /. Additionally, we assume 
that one particle is situated in £Ib close to the interface /. Assuming the particle 
crosses the interface in the first simulation step, a Dirac 5 function is created in fip 



close to the interface, as illustrated in Figure 3.4(a) This 5 function has a large 



impact on the region in £lp that is close to the interface, as the distribution p(x, At) 
is negligible in this region. Immediately after incorporating the particle into p(x, At), 
all its information is lost and we are forced to assume 100 independent particles with 
the probability distribution p(x, At) in Qp. In particular this implies that every 
particle has a 1% chance of being at the position of the S and 99% chance of being 
in the bulk distribution far away from the boundary. This is indeed not the case, 
since we know that there is exactly one molecule near the interface and 99 molecules 
away from the interface, but the nature of the continuum distribution means that this 
information must be lost or else we should necessarily demand a separate distribution 
for all molecules in flp. 

In the second diffusion step, we calculate p(x,2At) according to (3.6) and some 
'mass' a (2 At) may have drifted across the interface I (see Figure 3.4(b) ). Because the 
bulk distribution is far away from the boundary, almost all of this mass a(2At) conies 
from the 5 function close to /. Let us imagine that a new particle is now created in flp 
(according to the probability a(2At)), in which case the whole distribution needs to 
be rescaled, as shown in Figure 3.4(c) Because 99% of the mass in p(x, At) is situated 
in the bulk far away from the boundary, the majority of the rescaling happens in this 
region, such that effectively the mass needed to create the particle is almost entirely 
taken from the bulk, rather than from the region close to the interface. This also 
implies that most of the mass close to the boundary will stay and it is therefore 
possible to create another particle from this mass in further time steps. This is in 
contradiction to the result that would be expected if information was not lost in the 
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p(x,At) 

5{x-Xi{At)) 



'99\ 



1 
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p(x, At) 
p(x, 2 At) 



a(2At) n 



B 



X X 
(a) Particle jumps into Op and creates a Dirac (b) The distribution diffuses for one time step 

S function at the landing position. and an amount of the <5 flows back across the 

boundary. 



p(x, 2At) 
p(x,2At) 




>• Xi(2At) 



(c) If a particle is created, the majority of its 
mass is virtually taken from the bulk distribu- 
tion and the peak near the boundary remains. 

Fig. 3.4. Thought experiment that leads to errors in the variance. 



first time step. That is, the distribution close to the interface should dissapear and 
the bulk far from the interface is left alone. This effect is the main reason a higher 
than expected variance can be measured in fi^- Note, however, that this does not 
effect the expected values, as shown in Theorem |3.1| 

Of course, our thought example is an extreme case: we would expect that in prac- 
tical cases there would be a significant density of particles throughout the continuum 
region (otherwise we would be tracking them individually) . Nevertheless the fact that 
all information about an individual particle is lost as soon as it crosses the interface 
does generate an error in the variance of particle numbers near the interface, and the 
effect becomes more pronounced when the concentrations in Qp close to I are low. 

3.4.1. Dependence of the variance on the system parameters. We want 
to quantify the error in the variance as a function of the system parameters, which 
are the size of the domain [—L,L], the diffusion constant D, the simulated time T nna i 
and the total mass N. After a nondimensionalisation the macroscopic PDE (3.1) can 
be written in the form 



dn 
~dt = 

where the simulation is run until 



d 2 n 
dx 2 



x G [-1,1] 



1 final 



final 



D 



Hence, the system only has two parameters that need to be investigated: the final 
simulation time Tg nal and the total number of molecules N. 
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Fig. 3.5. Mean values and standard deviations of the number of iV/;('/j* •) depending on TZ , 
and N. Dashed line: expected number of particles; shaded area: expected standard deviation; dots: 
measured mean values; error bars: measured standard deviations. Other parameters are chosen as 
for Figure \3.3\ 



As in Figure |3^3| we use PBD algorithm (A.1)-(A.5) with n P = (-1,0), Q. B = 
(0, 1) and I = {0}. For each parameter (T£ Dai and N), we simulate the system 1000 
times for different values of this parameter and measure in each case the number of 



particles in fl B at the end of the simulation, i.e. the value N B (T£ aJ ). In Figure 3.5 
we see that the measured mean values match well with the expected outcomes. For the 
standard deviations, however, we see that for all values of Tg naI and N the measured 
outcomes are higher than expected. This is an undesired effect and the next section 
will discuss a way to reduce this artefact. 



4. A PBD algorithm with an overlap region. In Section 3.4 we saw that 
the immediate return of particles from fl P into Q B in combination with relatively 
low concentrations close to the interface can lead to errors in the variance of particle 
concentrations in fl B . One way to overcome this problem is the introduction of an 
overlap region where BD simulation and a continuum description exist in parallel, i.e. 
we will consider the case [B] defined in Section [2] by fijn flp ^ 0. A sketch of this 



new set up can be seen in Figure 4.1 with the overlap region denoted as O — f^nfip. 



We also denote the interfaces 1\ and I2 by 

h = dn B nn P , i 2 = n B ndn P . 

In the case of diffusion only, this new setup requires only subtle changes in the al- 
gorithm. Molecules are now incorporated into the concentration when they cross 



the boundary I\. The definition of p is still equal to (3.6 1, but we have to redefine 
a(t + At) and P2(x, t + At) as follows 



a(t + At)= p(x,t + At)dx, (4.1) 

Jn B \n P 

p 2 (x,t + At) = P{x > t +?f ) , brxen B \n P . (4.2) 
a(t + At) 

The introduction of the overlap region prevents undesired effects generated by molecu- 
les crossing over and coming back straight away, as discussed in the thought exper- 



iment in Section 3.4 In particular, a molecule initialized as a Dirac 5 function in 
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p(x,t) 


: N B (t) 






\Xi(t) Xj(t) 














t 




OfC Mb 





X 

Fig. 4.1. Sketch of the PBD algorithm with overlap region and the notation related to it. In f!p, 
molecules are described by their density distribution p(x , t) , in the microscopic domain Qb described 
by the number Ng(t) of molecules and their positions Xi(t), i = 1, . . . , Ng(t). In the overlap region 
O, both descriptions exist in parallel. The interfaces between the various subdomains are denoted I\ 
and I2 ■ 



(Bl) Calculate p(x,t + At) using (pDif and a(t + At) using (141) 



(B2) Generate uniformly distributed random number r in (0, 1). 

(i) If r < a(t + At), then create new particle in fl B \£lp according to 
the probability density P2(x,t + At) defined in (4.2). Set (3 — fiu\ 
where (3^ is given by (3.11). Set Nba = 1- 



(ii) If r > a(t + At), then set (5 = where is given by (3.11 ). 
Set N B ,i = 0. 

(B3) Compute positions Xi(t + At), i — 1,2, .. . ,Ns(t), of BD particles ac- 



cording to ( 2.4 1. 



(B4) Compute new concentration in flp by 

p(x,t+At) = /3p{x,t+At)+ S(x-Xi(t+At)) , for x G Op . 

xi(t+At)efip\n B 

(B5) Update the number of BD particles by 

N B (t+At) = N B (t)+N BA -\{x l (t + At) en P \n B ,i = l,...,N B (t)}\ . 

Terminate computation of trajectories of BD molecules which landed in 
Hp \ fl B (i.e. the BD particles which satisfy Xi(t + At) € f2p \S} B .) 
Then continue with step (Bl) for time t + At. 



Table 4.1 

One time step of the PBD algorithm with overlap region for system of diffusing molecules. 



Clp \ Cl B initially contributes very little to the overall probability density near the 
interface I2', by the time it has a significant probability of crossing I2 its distribution 
has become sufficiently spread that it is 'lost' in the subdomain VLp as is required for 
the continuous distribution. One time step of the second PBD algorithm is presented 
in Table |4.1| As before, we consider that the time step At is chosen so small that 
a(t + At) <C 1. Therefore, we only need to implement cases (i)-(ii) in step (B2), 
because the probability that two or more molecules are initiated in fl B during one 
time step is negligible. 
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Fig. 4.2. Simulation results of a diffusion process in Q = (— 1, 1) with no flux boundary con- 
ditions and initial conditions n(x,0) = 100<5(x + 0.95) with f2p = (—1,0), £Ib = (—0.1,1) and 
O = (—0.1,0) averaged over 1000 realisations. Dashed (red) line: expected outcome, (a) Solid line: 
mean value in Qp; bars: particle concentrations in Qb \ O (gray) and O (green), (b) (Gray) bars: 
measured variances in particle concentrations. Parameter values as for Figure 13. . 



In order to highlight the advantages of the overlap region, we simulate the same 
diffusion process as in Figure [3~3| with ftp = (—1,0) and Q B = (—0.1, 1). Then the 
overlap region is O = ftp fl ftp = (—0.1,0). The results are shown in Figure 4.2 



As before, the mean outcome matches well with the exact solution, with stochastic 



fluctuations inside the overlap region O due to the mixed description. In Figure 4.2(b) 



we see that the introduction of the overlap region indeed reduced the problem of high 



variances inside tip \ Op. A proof similar to Theorem 3.1 can be used to show that 
the PBD algorithm (B1)-(B5) satisfies Conditions (C.l) and (C.3) for At > 0. We 
will now further show that the algorithm describes a diffusion process exactly in the 
limit At -> 0. 

Theorem 4.1. Suppose that a PDE-description of the system is used in Qp = 
(—oo,0), and a BD simulation in flp = (—d,oo) with the overlap region O = (— d, 0) 
where d > 0. In the limit that At — > the expected concentration distribution 
Pp{x,t) = ¥\p{x,t)\ in the continuum regime and the expected concentration dis- 
tribution Pp(x,t) in the BD regime obey the equations 



dP P 
dt 


= D *Pr 
dx 2 


ox 


dP B 


= d &Pb 
dx 2 


- < p " 

ox 


dt 



5 (x + d) , x £ flp , 



5 (x) , x e ftp 



(4.3) 
(4.4) 



where 



P P {0,t) = 0, P B (-d,t) = 0, 



for t > 0. 



Extend each distribution to the whole line by defining Pp(x,t) = for x € (0,oo) 
and Pp(x,t) = for x € (-co, — d). Then the sum of these two processes n(x,i) = 
Pp(x,t) + Pp(x,t), x G Q, t > Q satisfies the diffusion equation (3.1 1. 
Proof. Consider the function 



n(x,t) = P P {x,t) + Pp{x,t). 



(4.5) 
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Fig. 4.3. Mean values and standard deviations of the number of particles in Qg \ f!p at 
time Tg nal depending on T^ nal and N. Dashed line: expected number of particles; shaded area: 
expected standard deviation; dots: measured mean values; error bars: measured standard deviations. 
Parameters as for Figure \4-S\ 



Clearly 



dn £)® 2n 
dt dx 2 



Since each of Pp and Pp is continuous at x 
continuous there. Moreover, since 



for x e (-00, -d) U {-d, 0) U (0, 00). 

— d and x = 0, function n will be 



dP F 



1 x=0 4 



dx 



dP F 
dx 



(0_,t), 



dn/dx is continuous at x — 0. Similarly, since 

-d+,t) , 



dPp 



dx 



Ox 



8P B 
dx 



dP B 
dx 



a;=0 4 



x=0- 



dP F 
dx 



(0-,t), 



=-d_ 



dP E 
dx 



{-d+,t), 



dn/dx is also continuous at x = —d. By standard regularity results this is enough to 
guarantee that n satisfies the diffusion equation on the whole real line. □ 
Remark. For the overlap region to give a different result from the simple PDB 
algorithm (Al)-(A5), we should choose d much bigger than the mean displacement 
\/2DAt of a particle given the time step At as defined in (2.4 1. 



4.1. Dependence of the variance on the system parameters. In order to 
show that the variances are indeed accurately produced by the PBD algorithm (Bl)- 
(B5), we repeated the numerical experiments conducted in Section 3.4.1 We chose 



ftp = (—1,0), O = (—0.1,0) and ftp = (—0.1, 1) and measure the number molecules 
situated in tip \ flp = [0, 1) at the end of the simulation. We again calculate the 
mean values and standard deviation and present the results in Figure 4.3 We can 



clearly see that the adjusted algorithm produces more accurate standard deviations 
than the original algorithm, whilst keeping the mean values correct. We conclude 
that the algorithm (B1)-(B5) produces an accurate BD simulation inside flp \ fip 
and therefore satisfies the Conditions (C.1)-(C.3). We will use the PBD algorithm 
(B1)-(B5) in the remainder of this paper. 
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5. Reaction-diffusion systems. In the next step we introduce chemical re- 
actions into the system presented in Section [4] We will concentrate on zero-order 
and first-order reactions |12) . First-order reactions are reactions which only have one 
reactant, for example, 

Xi X 2 , or X x A Xi + X 5l (5.1) 

where Xi denote chemical species and k\ (resp. k 2 ) the corresponding rate constant 
(which has physical units [sec" 1 ]). In what follows, we will denote by chemical 
species which are of no interest to a modeller. Then, considering that Xi is the only 



chemical species of interest, we can rewrite the reactions (5.1 1 as 

Xi 0, (5.2) 

where kd = k\ + ki . We will also consider zero-order reactions. An example is: 

-% X u (5.3) 

where the rate constant k p has physical units [Msec -1 ], i.e. it is the production rate 
per unit of volume and unit of time. It is relatively straightforward to implement 
zero-order and first-order chemical reactions in the PBD algorithms (Al)-(A5) and 
(B1)-(B5), because these reactions can be treated in the individual parts of the system 
(continuum and BD simulation) independently. Note that for higher-order reactions 
this is not necessarily the particles could react with the continuum inside the 

overlap region O. This case is not discussed in this paper. 

In the continuum regime reactions are represented by the term Rj(px,p2, ■ ■ ■ ,Pm) 
on the right-hand side of the reaction-diffusion P DE (|2.1[ ). For example, if the chemical 



species X\ is subject to chemical reactions (5.2 ) ( 5.3 1, then the reaction-diffusion PDE 



(2.1) takes the form 



dpi n d 2 p 1 
- W =D 1 —-k d p l +k p . 

In the BD simulations, the molecules act independently and the reactions can therefore 
be treated individually. A summary of how to implement various reactions in BD 
simulations can be found in |12) . 

Although implementation of the zero-order and first-order reactions is relatively 
straighforward, one has to still consider some special effects that are related to the 
coupling of the two parts of the domain. First, what happens when a particle that 
is supposed to react at time t\ crosses the interface I\ at an earlier time t% < t\t 
Since we assumed that all information about particles is lost as soon as they cross the 
interface I\ , we incorporate it into the continuum and the reaction at time t\ does not 
happen. Second, what happens when a particle is created inside the overlap region 
Ol A number of solutions to this problem are possible: one could split the creation 
in equal parts, or declare creation to only contribute to either the continuum or the 
molecular-based description. We will here assume that all creation inside the overlap 
region occurs in the form of molecules with exact positions. 

Finally, let us note that (reactive) boundary conditions on the external boundary 
<9£1 can be treated according to the corresponding modelling regime, i.e. whether the 
corresponding segment of dfl is part of dtlp or dfts- Derivation of reactive boundary 
conditions of BD simulations which are consistent with the PDE description can be 
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found in [8!. External boundaries slightly modify the computation of p{x,t + At) in 
Op. It is still given by (3.6) but the kernel (3.7) has to be updated to take into account 
the boundary condition imposed on the external boundary 90. We have already done 



this when we showed simulations of the diffusion process in Figures 3.3 and 4.2 in the 
hnite interval (—1,1) with no flux boundary conditions. However, for small timesteps 
At the change is negligible, since all the action takes place near the interface 1-2.. 

We conclude this section with three examples which illustrate the behaviour of the 
PBD algorithm (B1)-(B5) for reaction-diffusion systems. They include the modelling 
of morphogen gradients and chemisorption. 

5.1. Example 1: morphogen gradient. In the first example we compute a 
steady state for a morphogen gradient model [28], [32| [20] . We consider one chemical 
species (morphogen) inside the domain O = ( — 1, 1). All parameters are dimension- 
less for simplicity. The only reaction inside O = (—1,1) is the degradation (5.2). 



Additionally to this reaction, we assume a constant influx J/D through the left-hand 
boundary x = — 1 to the continuum subdomain Op = (—1,0). We use O = (—0.1,0) 
and Op = ( — 0.1, 1) with a no flux boundary at x = 1. Since we only have first-order 



reaction (5.2), the exact solution n{x,t) is given by 



dn 

~dt 



D 



d 2 n 
dx 2 



kin , 



dn 
dx 



(-M) = -J, 



dn 

dx 



(M) = 0. 



(5.4) 



This system is initialised with n(x,0) = 0, x € O, and we let it run until it (approx- 
imately) reaches the steady state. We use J = 1000 and ki = 1. The second PBD 
algorithm (B1)-(B5) is run with the same parameters presented in Section [4] for time 
Tfinai = 20 which is (approximately) a time at which the model settles into the steady 
state. The reaction (5.2) is simulated in a time-driven manner in Op, which means 



that for each morphogen molecule it is decided randomly at the end of each time step 
whether it was degraded or not (the probability of degradation of each molecule is 
equal to kd At provided that kd At <C 1). 

The result of a single simulation of the PBD algorithm (B1)-(B5) is plotted in 
Figure 5.1(a) We plot the PDE solution in Op \ O as a black line and the (gray) 
histogram of molecules in Op \ O. In the overlap region O, we compute the total 
"number" of particles by 



N (Tt 



final 



= / p(x,T t 
Jo 



fi„ai)dx + I {xi(T fina] ) 60,1 = 1,2,..., A^(T final )} I. (5.5) 



The value of No(Ta na x) is plotted as the green bar in Figure 5.1(a) The results of 



a single simulation of the PBD algorithm (B1)-(B5) are compared with the exact 
solution, which is obtained by solving the PDE (5.4) numerically until time Tg na i. In 



this simple example, it is also possible to find an analytical expression for the steady 
state profile which is approximately equal to the presented dashed line. 

Note that the jagged appearance of the continuum solution close to the interface is 
not numerical error, but represents the fact that as molecules cross from the discrete to 
the continuum side information about their exact location is lost gradually over time 
(remember that Figure 5.1(a) shows just one realisation of the stochastic process). The 
corresponding distribution on the discrete side Op \ Op would be 8 function spikes at 
the location of the particles, which we have in effect locally averaged by our binning 
process. Thus the jaggedness can be seen as a gradual transition in the solution from 
isolated discrete particles to a continuum density distribution. An ensemble average 
over 100 simulations of the PBD algorithm (B1)-(B5) is shown in Figure 5.1(b)| The 
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Fig. 5.1. Simulation results for Example 1. Dashed (red) line: exact solution given by ( |5.4| ; 
solid line: p{x,Tf lna x); (gray) bars: spatial concentration of particles at t = Tg na i in fig \0; (green) 
bar: No(Tfi n!i i) given by |5.5[ l. Parameters as described in the text. 



stochastic fluctuations are reduced compared to the single simulation and the results 



compare well with the exact solution for the expected probability density (5.4 1 



5.2. Example 2: reversed morphogen gradient. In this example we intro- 



duce a second reaction in addition to (5.2) - a local production of molecules 



.4 



1], 



(5.6) 



where x s defines the size of the creation zone. As before we consider all parameters 
to be dimensionless: k p is defined as the rate of production per unit length. For 
this system we use no flux boundary conditions on both ends. The combination 



of localized production (5.6) and degradation (5.2) ensures that the system settles 



into a non-trivial steady-state which we will compute with the PBD algorithm (Bl)- 
(B5). The exact solution (which the PBD algorithm (B1)-(B5) approximates) can be 
described by 



dn 



= D 



dx 2 



fan + faX[x„i] 



dn 
dx 



(-M) = o, 



dn 
dx 



(M) = o, 



(5.7) 



where X[x 3 ,i] ^ s ^h e characteristic function for the interval [x s , 1] that takes the value 
1 inside and outside of the interval. The production reaction (5.6) was implemented 



in BD simulations in an event-driven way, such that particles can get created at any 
time in-between time steps and the number of particles created in one time step is 
not limited. We used k p = 1, fed = 2000 and x s = 0.5. For the PBD simulations 
we use the same parameters as in Section [4] In particular, we have £lp = (—1,0), 
O = (-0.1, 0) and fl B = (-0.1, 1). 



A single realisation of this process is plotted in Figure 5.2(a) We plot the PDE 
solution in fip \ O as a black line. The concentration of molecules in £Ib \ O is 
visualized as a (gray) histogram. In the overlap region O, we plot Ao(Tfi na i) given by 
( 5.5 ) . The concentration gradient is now reversed, as the creation of particles happens 



near the right-hand boundary. Again, one can clearly see the stochastic fluctuations 
in this plot which also have effects on the value of p(x, t) far from the overlap region 



O. However, as we draw an ensemble average over 100 simulations in Figure 5.2(b) 
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Fig. 5.2. Simulation results for Example 2. Dashed (red) line: exact solution given by \5.7) ; 
solid line: p(^,7fl na i) in Qp \ O; (gray) bars: spatial concentration of particles at t = Tg na i in 
Qb \ O; (green) bar: No(Tfi n!i i) given by l |5.5| . Parameters as described in the text. 



the results converge towards the exact solution which is obtained by solving the PDE 



(5.7). It is plotted as a (red) dashed line in Figure 5.2 



5.3. Example 3: chemisorption. Our last example is the polymer coating of a 
virus surface [T31 E] • We will describe it as irreversible adsorption (chemisorption) of 
polymers to a two-dimensional surface as was introduced in [9| . This example presents 
a typical application area of PBD algorithms. A detailed model is used close to the 
reactive boundary where positions of individual molecules influence the dynamics of 
diffusion-driven adsorption. On the other hand, a less detailed model can be used far 
away from the adsorbing surface. In the bulk the behaviour of reactive polymers can 



be described by the macroscopic reaction-diffusion PDE (2.1| in the form 



dp d 2 p 

i = D dS~ kdP - (5 - 8) 

in the semi-infinite domain Q = (0, oo) (here, z is the distance from the reactive 
surface which is at z = 0). This equation takes into account two processes which 
mainly influence chemisorption dynamics 9 : diffusion of polymer molecules and the 
hydrolysis of reactive groups in the solution. Both processes can be implemented in 
the BD context as we saw in the previous examples. However, this level of detail is 
only needed close to the virus surface. 

Whenever a polymer molecule interacts with the surface, it is either reflected 
or (irreversibly) adsorbed. The chemisorption is modelled by a random sequential 
adsorption (RSA) algorithm [7]: we check whether the corresponding binding site 
on the surface is free and then the reaction occurs with a certain probability. This 
probability is related to the reaction rate constant of the binding reaction as given in 
[S]. The reader can find more details about the model in 9 . In this paper, we show 
that the PBD algorithms can be used to compute the results from [5] . We will use the 
same parameters, namely D = 5 x 10~ 5 mm 2 s _1 , k c i = 1.3 x 10~ 4 s _1 and At = 0.01 s. 



Then the mean displacement per time step accor ding to K4h is V2DAt = 10" 3 mm 
and we therefore choose the size of the overlap region of the PBD algorithm (B1)-(B5) 
as \0\ = 10 -2 mm. From the results in [5], we estimate that a maximum length of 
L = 2 mm is enough to simulate the binding process and use the Dirichlet boundary 
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Fig. 5.3. Example 3 (chemisorption to virus surface). (Gray) histograms: concentration profile 
in molecules/mm at two given times computed by the PBD algorithm (B1)-(B5); solid line: results 
of the RSA-PDE model presented in Parameters as shown in the text. 

condition 

n(L,t) = c Q cxp(-k d t) , 

where cq — 1.2 x 10 4 molecules/mm is the initial concentration of molecules (i.e. 
n(z, 0) = Co for z E £1). To apply the PBD algorithm (B1)-(B5), we choose fls = 
[0,1.01) mm, O = (1,1.01) mm and ftp — (1,2) mm. The RSA algorithm was per- 
formed using a nearest neighbour exclusion on a 100 x 100 grid of receptor binding 
positions on the surface [TJ. 

In Figure [573] we plot the concentration profile inside Qb of a single simulation at 
two different times (gray histograms) . We compare the results of the PBD algorithm 
(B1)-(B5) with the results of the RSA-PDE model presented in [5] (black lines). As 
shown in [5], the RSA-PDE model also compares well with the full BD simulation. 
The number of molecules which are attached to the surface as a function of time is 
plotted in Figure [5~4] (six realisations of the PBD algorithm (B1)-(B5) are plotted as 
green solid lines). Again, we see an excellent agreement with the result from 9 which 
is plotted as the black dashed line. 

6. Discussion. In this paper we have presented two PBD algorithms that com- 
bine Brownian dynamics with mean-field reaction-diffusion PDEs. This method pro- 
duces exact Brownian dynamics simulations in one part of the domain and couples 
them with mean field approximations in another part of the domain. An algorithm of 
this type is useful for various application areas in computational biology and beyond, 
for example, when a detailed description of individual molecules is required near a 
receptor or ion channel, but becomes impractical in the bulk of a cell [SJ 116) ; or when 
a detailed stochastic simulation of actin dynamics is required inside filopodia, but be- 
comes impractical in the bulk of a cell |33j . Another application area, chemisorption, 
was discussed in Section |5.3| By using our approach it would also be possible to use 
finite-sized particles in the BD simulation and couple these with the corresponding 
mean field results presented in J3J. 

In the literature several hybrid models have been developed in the context of fluid 
dynamics, but they do not discuss issues that arise from the incorporation of chemical 
reactions. Alexander et al pQ presents a hybrid model that uses virtual particles at 
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Fig. 5.4. Example 3 (chemisorption to virus surface). Number of polymer molecules which are 
bound to the virus surface as a function of time. (Green) solid lines: six realisations computed by 
the PBD algorithm (B1)-(B5); (black) dashed line: results of the RSA-PDE model presented in (Pp. 
Parameters as shown in the text. 



the grid point nearest to the interface to calculate fluxes across the boundary and 
to generate accurate density fluctuations inside the particle region. Reference |31j 
extends this approach by the introduction of an overlap region similar to O introduced 
in Section [4j An identical flux exchange with particles confined to a grid is presented 
in [17] . Chemical reactions in the solution were considered in [5]. This model was 



discussed in Section 5.3 It couples Brownian dynamics of molecules in the solution 
with a more detailed description of the adsorbing boundary. In [9] a hybrid (RSA- 
PDE) model has been developed which replaces BD in the solution by solving the 



PDE (5.8) with a suitable stochastic boundary condition. The PBD algorithms are 
able to replace the stochastic boundary condition by a (small) BD region close to the 
surface. Although the hybrid RSA-PDE model introduced in [3] was sufficient in the 
case of (irreversible) adsorption, the situation is becoming more challenging whenever 
the binding reaction is reversible |21j . In this case, a molecule which is released by 
the surface will initially stay close to the surface and can rebind to the same receptor 
(binding site). This geminate recombination can be captured by the PBD approach. 
Reversible reactions are common in biological applications [211 116] . 

Hybrid approaches for reaction-diffusion processes which couple different mod- 
elling approaches have also been introduced in the literature [23J |T5l [6] [13]. A 
mesoscopic lattice-based description coupled with macroscopic Fisher-Kolmogorov- 
Petrovsky-Piscounov PDE was used in [53] to study front propagation in a lattice- 
based reaction-diffusion model. A hybrid model for reaction-diffusion systems in 
porous media that combines pore-scale models with Darcy-scale models is presented 
in [57]. Flegg et al [TS] introduced the so-called Two Regime Method which couples 
a lattice-based (compartment-based) reaction-diffusion model with BD simulations. 
One advantage of the PBD algorithms over the Two Regime Method is that there 
are more efficient tools for PDE simulations than for compartment-based reaction- 
diffusion models. On the other hand, compartment-based models provide more de- 
tails (including fluctuations) and hybrid models which couple BD simulations with 
compartment-based models do not require the overlap region [HI [13] . Since it is 
possible to couple the macroscopic PDE description with (mesoscopic) compartment- 
based models [ja) and compartment-based models with (microscopic) BD simulations 
[15j . then an alternative approach to PBD algorithms would be to use compartment- 
based models in the overlap region. That is, the computational domain would be 
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divided into three regions where the PDE, compartment-based and BD descriptions 
would be used. These three regimes would be coupled using the results from the 
literature [51 [T5]. Compartment-based models and macroscopic PDEs can also be 
coupled through another intermediate regime using a tau-leap method |13j . In this 
paper, we showed that PDE models and BD simulations can be coupled without using 
intermediate compartment-based models. 

The PBD algorithms presented in this paper should be seen as a first step to- 
wards a more general setting. Some parts of the algorithm (B1)-(B5) extend easily 
(at least theoretically) into higher dimensions, but in practice additionally difficulties 
are posed. One example is the necessity to sample from a multidimensional proba- 
bility distribution to find the position of newly created molecules. Additionally, in 
higher dimensions one can also expect to deal with higher order reactions, including 
bimolecular reactions. For a discussion of how to implement bimolecular reactions 
for BD simulations, we refer to [TO], but the real problem occurs inside the overlap 
region O = Qb H flp, where a molecule could react with another molecule, or with 
the continuum. Since the reaction-diffusion PDEs are solved numerically using a suit- 
able mesh, it is important to study methods for coupling individual molecules with 
the numerical discretization of macroscopic PDEs [18]. For the ion channel applica- 
tion mentioned before, one also needs to think about how to incorporate electrical 
charges and resulting forces into the system. One can imagine that these forces act 
as a boundary condition on the continuum model and as an effective force on the 
particles. 
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